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ABSTRACT 

We suggest a novel discretisation of the momentum equation for Smoothed Particle 
Hydrodynamics (SPH) and show that it significantly improves the accuracy of the 
obtained solutions. Our new formulation which we refer to as relative pressure SPH, 
rpSPH, evaluates the pressure force in respect to the local pressure. It respects New- 
tons first law of motion and applies forces to particles only when there is a net force 
acting upon them. This is in contrast to standard SPH which explicitly uses Newtons 
third law of motion continuously applying equal but opposite forces between particles. 
rpSPH does not show the unphysical particle noise, the clumping or banding instabil- 
ity, unphysical surface tension, and unphysical scattering of different mass particles 
found for standard SPH. At the same time it uses fewer computational operations, 
and only changes a single line in existing SPH codes. We demonstrate its performance 
on isobaric uniform density distributions, uniform density shearing flows, the Kelvin- 
Helmholtz and Rayleigh-Taylor instabilities, the Sod shock tube, the Sedov-Taylor 
blast wave and a cosmological integration of the Santa Barbara galaxy cluster forma- 
tion test. rpSPH is an improvement these cases. The improvements come at the cost 
of giving up exact momentum conservation of the scheme. Consequently one can also 
obtain unphysical solutions particularly at low resolutions. 

Key words: Numerical Methods-Smoothed Particle Hydrodynamics- 
Hydrodynamics-Instabilities 



1 MOTIVATION 

The smoothed particle hydrodynamics method was invented 
by iLucyl (119771) and iGingold and Monaghanl (119771), both 
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with interests in astrophysical applications. Besides an enor- 
mous literature of successful application also many short- 



comings of it have been presented in the literature 1 


Stein- 


metz and Mueller |1993| |Herant 1994| Swegle et al. 


1995| 


Imaeda and Inutsuka 20(52) Agertz et al. 2007; Read et al. 


2009 1 often suggesting a fix to the reported problem 


Mon- 


aghan| |1994| |Cummins and Rudman 


|1999| |Rasio| 


2000 


Attwood et ah] 2007 Hu and Adams 


2007[ |Graham and 


Hughes|2008| Price|2008j B0rve et al.|2009||Rafiee and Thi- 


agarajan|2009| Xu et al.|2009| to name but a few) Similarly 



there have been troubling news of how seemingly small dif- 
ferences in the initial setup led to very unexpected results 
(e.g. |Lombardi et al.]|1999[ ). This is all somewhat surpris- 
ing given that in many of the cases where large inaccuracies 
have been found the only relevant equation (besides mov- 
ing the particles f = v) stems from the pressure gradient 
accelerations 
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where D/Dt denotes the Lagrangian derivative, and p the 
pressure. In what follows we describe a new discretization 
of the momentum equation that avoids essentially all of the 
previously known problems of SPH. We will refer to this 
new method as "relative pressure SPH" or abbreviated as 
rpSPH. We will first describe it, discuss implementation de- 
tails and then present results for relevant tests highlighting 
the superior performance of this new approach. 

All the simulations shown here are carried out with 
Gadget-2 (Springcl 2005) (version 2.0.4) with only most mi- 
nor changes explained in the text. The appendix describes 
how to convert Gadget to our rpSPH formalism. 



2 RELATIVE PRESSURE SPH 



The equation of motion without viscous or gravitational 
forces in essentially all SPH codes and Gadget-2 ( |SpringeI| 
and Hernquist|2002 Springel 2005), in particular, is 
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where the fi are defined by 



fi 



1 + 



hj dpi 
3pi dhi 



(2) 



(3) 



and the abbreviation Wij(h) = W(\fi — fj\, h) has been used 
for the kernel function, W. It employs variable smoothing 
lengths so that the number of neighbours for each particle 



with 



^ hi is maintained at a nearly fixed value iV S ph. 



The compact cubic spline kernel is used which is summarized 



Monaghan (19921, extends to radii as large as the smooth- 



ing length h and is zero outside. While many choices would 
exist to use a different discretisation here flMonagha n|1992 



Rosswog 2009) most previous work we have found essen- 
tially retains a form very close to equation |2| or uses an- 
other symmetric form that sums over (Pi + Pj)(piPj). The 
reason previously given for these choices are their symmetric 
form encapsulating Newtons third law of motion, the action- 
reaction law. By guaranteeing that particles give pairwise 
identical but reversed forces one ensures linear momentum 
conservation of the entire scheme. The key here is that par- 
ticles are always pushing as soon as they have any pressure 
regardless of whether there is a pressure gradient. The ones 
with the highest pressure values are pushing the most. When 
one has a large number of particles in a perfectly symmetric 
configuration all the pushing will average out for an indi- 
vidual particle. This is to some extent what happens in real 
gases. The pressure itself is mediated by the collision of the 
molecules the gas is made of. 

From a physical point of view of a Lagrangian fluid 
element, however, one should only be interested in the pres- 
sure forces of neighbouring fluid elements exerted on oneself, 
since the actual equation of motion is pv — VP. This is the 
main idea of rpSPH, a particle is accelerated only if a force 
is acting upon it, i.e. Newtons first law of motion. rpSPH 
derives its equation of motion directly from equation |2| 
by subtracting the constant pressure of the particle under 
consideration from the pressures of all the particles being 
summed over. Since a gradient is computed, the subtraction 
of a constant does not change the mathematical meaning of 
the difference equation. However, as we will demonstrate it 
dramatically affects the error properties of the entire scheme. 
The resulting equation of motion reads 



dvi 
df 



=-£< 



3=1 



P, 



V l W ij {h i ) 



(4) 



One immediately notices that this formulation breaks the 
symmetry between the pairwise forces of particles. Particles 
that have a pressure difference are both accelerated into the 
same direction along the pressure gradient. Linear momen- 
tum conservation hence will only be achieved if the modeled 
pressure gradients are resolved. On the other hand if one 
does not resolved the relevant pressure gradients one cannot 
possibly get a correct solution to a hydrodynamic problem 
in any case. 

After all, it is important to recall that when construct- 
ing conservative schemes one does not necessarily minimize 
the numerical errors but rather ensures that one is making 
symmetric errors so that the conserved quantity does not 



change. Consequently, in rpSPH monitoring the total angu- 
lar and linear momentum is an indicator of whether one may 
have resolved the relevant length scales. 

Many of the advantages of the entropic function based 
SPH formalism Spring el and Hernquist| ( |2002[ > stem from 
avoiding the PdV term that generally is discretized analo- 
gous to equation B. So in this formalism rpSPH is partic- 
ularly trivial to implement. It involves setting the first term 
on the right hand side of equation |2| to zero and change 
the second by subtracting the pressure of the particle under 
consideration. This literally is achieved by modifying one 



line of code in Gadget-2 Springel ( 2005 ) as shown in the last 



appendix. The resulting scheme saves two multiplies, one 
division and one addition for one additional subtraction in 
the main loop over neighbors. So there is no performance 
penalty in using rpSPH as compared to standard SPH. 

rpSPH is seemingly close to equation (3.1) of |Monaghan| 
(1992) first discussed by |Morris| ( |1996a| ) which we will refer 
to as the Morris formulation. It reads, 



dvi 
diT 



N 
3=1 



Pip j 



(5) 



Monaghan dismissed his version for two reasons. The 
first is that "an isolated pair of particles with different pres- 
sures would bootstrap themselves to infinity" and the sec- 
ond is that it is difficult to construct a consistent energy 
equation. The latter is irrelevant in the formalism evolving 
an entropic function |Lucy| ( |1977[ ); Springe !" and Hernquist| 
( 2002 1 in which the PdV work does not enter. The first 



reason we find unappealing since it is actually the correct 
solution. The simulation having two particles estimates a 
pressure gradient. So over the model volume, i.e. the two 
particles and their smoothing volumes there exists a mono- 
tonic pressure gradient. Both particles hence should be ac- 
celerated along it. Interestingly, Monaghan did not discuss 
the equivalent case for the symmetric standard SPH. In this 
case both particles push each other to infinity no matter 
what. If they have the identical initial pressures their cen- 
ter of mass will not change if they vary their center of mass 
moves exactly as in rpSPH. In rpSPH they will move to- 
gether while in SPH they will accelerate each other apart to 
infinity. We have tested this on a spherical blob of material 
in vacuum. We set the pressure of the particles after the 
densities have been computed from kernel smoothing. This 
way all particles have identical initial pressure. The config- 
uration is completely stable in rpSPH yet blows itself apart 
in SPH in just a sound crossing time. The reason why we do 



not choose equation (3.1) of |Monaghan ( 1992 \ is because we 
find it to be unstable at least with the leapfrog time integra- 
tor in Gadget-2 (see Figure [5] below). Another formulation 
close to rpSPH discretization we could find in the literature 
is presented Morris et al.| ( |1997a[ ), who chose to subtract a 
background pressure. This still leaves an equation of motion 
in which the pressure of the particle under consideration 
remains part of the hydro force estimate. 

An easy way to see why our discretization is valid (Wad- 
sley, 2010, private communication) recognizes rrij/pj = dV 
as the volume element dV and sees that equation [4] is 
equivalent to j [VPp -1 — PVp -1 ] dV , which is the same 
as J p~ VP and is the term we want. The previous version 
discussed by Morris | ( |1996a l in contrast is the discretized 
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form p _1 VP rather than our VPp -1 - PVp -1 . Note that 
our form is also not the general form suggested by equa- 



tion 2.13 of Monaghan (20051 and in this regard is a new 



formulation. The striking aspect is that in the actual differ- 
ence form all that is new is that one index that used to be 
i is now j. So this literally is a one letter change to codes 
that implement the Morris formulation. How this can lead 
to a dramatic change in accuracy becomes obvious from the 



standard error analysis. Price (20041 gives the error of the 



summation interpolant in equation (3.9) and the error of the 
gradient opera tor in equation (3.11). In a formulation as by 



Morris 



( 1996a I one discretises p i *VP and sees that the er- 
rors in the interpolation of p _1 and the errors of the chosen 



VP discretisation multiply. What Morris ( 1996a ) realized 



was that there are no error terms in his discretisation for 
constant functions in the gradient operator. However, the 
complete error terms still end up being the product of the 
density estimate and the pressure gradient estimate. The 
advantage of the rpSPH discretisation is that its error terms 
are the one of a gradient and are not further multiplied by 
errors of a density interpolation. It also retains the vanishing 
error for constant functions of the Morris discretisation. 
Interestingly, a linear stability analysis reveals that rp- 



SPH has the same dispersion relation as the form of ( Morris 



1996a his equation 10). He has shown that this form 



• is always stable, independent of the background pres- 
sure, 

• has a numerical sound speed that depends less on the 
particle spacing as compared to standard SPH, 

• and does not have unphysical unstable transverse waves 
in two nor three dimensions when using kernels with com- 
pact support. 

rpSPH retains all of these advantages while at the same time 
reducing the discretisation error. 

In the standard SPH approach there are in fact in- 
finitely many possible choices for the discretisation of the 
pressure equation (equation 3.5 in Monagh an|1992 1. This is 
also true for rpSPH. E.g. P^V P CT /p- PVl/p = ap _1 VP 
which suggest the discretisation 



dvi 
~dt 



N 

E 1 

3=1 



(6) 



for any a different from zero. We have verified that many 
choices of a work for a variety of test problems. In the follow- 
ing, however, we restrict our attention to the case of a = 1. 

Whether these theoretically advantageous properties of 
rpSPH hold up in practice is assessed in a range of test 
problems in the following section. 



3 TESTS OF RPSPH 

We will employ a Courant factor of 0.3 (i.e. 0.15 in Gadget 
where the kernel has a maximal radius of h). 



3.1 Reduced Velocity Noise 

We start with 50 2 particles on a periodic regular lattice with 
7 = 1.4, a sound speed and uniform density of unity and zero 
initial velocities. The particles should stay at rest. However, 



as we can see in Figure [T] the total kinetic energy in the vol- 
ume grows rapidly. The total energy in the system is, how- 
ever, conserved to better than 5 x 10~ 5 of the initial value 
for these tests at a value (7(7 — 1)) _1 ~ 1.78. So the kinetic 
energy growth in the particle distribution only corresponds 
to about less than one in one thousand of the total. I.e. the 
kinetic energy the particles obtain is taken from a slightly 
decreasing internal energy allowing the total to be conserved 
to high precision. The lower the neighbour number the faster 
that growth. The maximum noise reached is controlled by 
the artificial viscosity. The noise also decreases only very 
slowly over time after reaching the maximum. This is one 
of the main reasons why particle settling is so important in 
SPH simulations. The slow decline also shows why in general 
settling procedures can be computationally quite intensive. 
In the same figure we also plot results using rpSPH which 
dramatically reduces this spurious kinetic energy creation 
keeping it at zero to machine precision. Rasio (2000) caution 
that it makes no sense in standard SPH to increase parti- 
cle numbers while keeping the number of neighbors fixed. 
Once a neighbour number is reached that keeps noise in the 
force calculation to a minimum we find rpSPH to be stable 
while only increasing the particle number. Note that we also 
have ran these tests dramatically reducing the Courant fac- 
tor without any improvement in the case of standard SPH. 

The thick solid line in Figure [I] uses 20 neighbors which 
seems optimal for this 2D calculation with the cubic spline 
kernel. Here one has enough neighbors to estimate the gra- 
dients more accurately while still having too few neighbours 
to show its pairing instability. So one may be tempted to dis- 
miss the finding that one has the large velocity noise as long 
as one uses the "correct" number of neighbours in one simu- 
lation. Unfortunately, this best choice, however, is only ap- 
plicable at the uniform density. To show this we perturb the 
x positions by a small amount so that the initially uniform 
Xq positions are changed by adding sin(2-7r a;o)/25 to them 
which gives central densities that are about 30% above the 
mean. We keep again the pressure to be exactly constant by 
setting the entropy of the gas only once the density has been 
estimated from kernel smoothing. The thick long dashed line 
in Figure [T] gives the associated velocity noise. It again is of 
order one percent of the sound speed and grew very rapidly. 

One may also be tempted to to dismiss this particle 
noise as irrelevant as it only contains less than a tenth of 
a percent of the total energy of the system. However, we 
will see in the following it is what leads to unphysical shear 
viscosity once one considers shear flows further below. 

3.2 Absence of the Pairing Instability 

The velocity noise we just discussed is unfortunately not 
isotropic nor is it random. It has a dominant component 
for velocities towards directions of other particles and is an 
effect that aids the pairing instability. Our reasoning here is 
somewhat contrary to the explanations in the literature as 
e.g. in 



Schuessler and Schmitt| ( |l98l[ ), |Vaughan et al.| ( |2008 l 
Read et al. ( 2009[ ). These studies suggest that it is the 



shape of the kernel function that causes clumping instability. 

Given that rpSPH does not show any spurious velocities 
in the uniform density test given above while standard SPH 
develops clumps within a few sound crossing times already 
shows that it cannot be the shape of the kernel alone that 
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Figure 1. The total kinetic energy in the static uniform density 
test as a function of time. The labels give the square root of 
the number of particles used, the number of neighbors and the 
artificial viscosity parameter. So the 50-30-aO.l used 50 2 particles, 
with 30 neighbors and and artificial viscosity parameter a = 0.1. 
We clearly see that using more neighbors leads to a slower growth 
of the unphysical kinetic energy. Using higher viscosity values 
can not stop that growth but does influence how quickly this 
noise is damped over time. rpSPH even in the lowest numbers 
of neighbors reduces the kinetic energies associated with these 
unphysical motions by at least seven orders of magnitude and for 
sufficient neighbor numbers keeps it at zero to machine precision. 



is relevant here. In the following tests we have looked for 
any sign of the clumping instability but have not found any 
evidence for it independent of the number of neighbors we 
used. This fact is certainly one of the factors that makes our 
formulation have radically reduced errors in general. 

The clumping instability stems from particles pushing 
each other closer to other particles. With a smaller distance 
to the other particles the gradient of the kernel becomes 
smaller and in the next time step the particle gets pushed 
further away from its initial position. This way particles can 
pile up in the flat central part of the smoothing kernel. That 
rpSPH does not show the clumping instability makes one 
hopeful that even higher order kernels could now be used to 
further improve the accuracy of the obtained solution. 



3.3 Dramatically Reduced Numerical Shear 
Viscosity 

An easy two dimensional setup uses an adiabatic index of 
7 = 1.4 in a unit square domain x € {0, 1}, y G {0, 1} 
with periodic boundary conditions. Particles are initialised 
on an exactly square lattice with a density of p(x, y) = 1 so 
that the initial density estimate from the SPH kernel gives 
a density estimate of unity to better than four parts in one 
thousand. We then add a sinusoidal velocity perturbations 
to this uniform distribution. We set the pressure to Po = p/j 
to have a sound speed of unity. For the first tests here we 
only use 50 2 particles as there are no features to resolve. In 
all cases we evolve to time t — 4. 



shear flow setup gives an average kinetic energy of 1/8. A 
detailed discussion of how SPH behaves on this test for dif- 
ferent neighbor numbers, particle numbers and viscosity pre- 
scriptions is given in the Appendix. In summary, it does very 
poorly and transfers kinetic energy into heat very rapidly 
loosing tens of percent in only four sound crossing times 
(two crossing times of the fastest particles). It also gives 
more dissipation when using more particles which makes a 
convergence study at fixed neighbour number impossible. 
Below, when we discuss rpSPH for viscous flow, we show 
that the effective numerical viscosity of standard SPH is 
non-Newtonian and very large which explains why standard 
SPH is inadequate to modeling fluids in general. 

The results for rpSPH are summarized in Figure[2]which 
is to be compared to the bottom panel for SPH, plotting 
the kinetic energy in the box as a function of sound crossing 
times. Note that the y-axis in the two panels of Figure [2] is 
different by a factor of 30. 

We should note that this test problem when used in 
typical Cartesian grid codes will give zero numerical dissi- 
pation to machine precision because the uniform nature of 
the flow along the grid axes. 



Price (20041 considered seemingly similar tests in his 



thesis. However, note that there an isothermal equation of 
state was used. He shows one case with c s = 0, i.e. a pressure 
less fluid and another case with c„ — 0.05. His initial shear 
profile is like the one we choose here but with twice the am- 
plitude so the fastest particles move with unit speed. The 
pressure-less case is irrelevant here since without pressure 
forces the discretization of the momentum equation cannot 
make a difference. For the second case with c s = 0.05 he 
gives the result only after one sound crossing time at t = 20. 
The density fluctuations have grown to order of a percent 
of the initial value after that single sound crossing time. As 
Figure [2] shows for most reasonable choices of neighbor num- 
bers not much kinetic energy is dissipated over this time also 
in our simulations with an adiabatic index of 1.4. We also 
repeated this isothermal shear test and find results consis- 
tent with | Price | ( |2004[ ) . This again emphasises that as long 
one is interested in very few sound crossing times, SPH can 
give correct results and gives an indication that this is not 
specific to Gadget. 

This problem also allows us to measure the effective 
Reynolds numbers one can hope to model with standard 
SPH. 

Solving analytically the incompressible Navier-Stokes 
equations for our initial conditions because to good approx- 
imation only the viscous term is relevant we have 

dv 2 

Since the variables are separated we can easily find that 



v(x, t) — e 



' v(x, t = 0), 



(7) 



showing that the initial shape of the profiles should not 
change. We can now use equation [7] to get an estimate of 
the kinematic viscosity from the fraction of kinetic energy 
remaining 

Ul/F) 
8n 2 t ' 



(8) 



We choose v x (y) = Sv y cos(2ny) with Sv y — 1/2. This where F denotes the fraction of the kinetic energy remaining 
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Figure 2. The fraction of the total kinetic energy in the uniform 
density shear test as compared to its initial value as a function of 
time for rpSPH. The labels give the square root of the number of 
particles used, the number of neighbors and the artificial viscosity 
parameter. So the 50-40-aO.l used 50 2 particles, with 40 neigh- 
bors and an artificial viscosity parameter of a = 0.1. We clearly 
see that using more neighbors always leads to less artificial shear 
viscosity. With more particles the effective numerical dissipation 
(which is very small to begin with) decreases. The corresponding 
panel for SPH below has a j/-axis 30 times as big. (Bottom) The 
fraction of the total kinetic energy in the uniform density shear 
test as compared to its initial value as a function of time for stan- 
dard SPH. We clearly see that using more neighbors always leads 
to less artificial shear viscosity. However, increasing the particle 
number leads to more dissipation. 



up to time, t. This is v ~ 3.2 x 10 -3 ln(l/F) at four crossing 
times. 

The Reynolds number measures the ratio of inertial 
forces, pv 2 , to viscous ones, fj,V/L, where L is the charac- 
teristic length scale, V the mean velocity and p is the dy- 
namic viscosity. So R — VL/v with v — p/ p the kinematic 
viscosity. Despite the ambiguities we may take L — 1/4 



the quarter wavelength of the velocity perturbation, V = 
JNp 1 J2 N v 2 = 1/2^2 the root mean squared velocity. 

So for the typical values of F we found for SPH say 
70 (97) per cent remaining equation |8| gives a kinematic 
viscosity of ~ 10~ 3 (1CP 4 ). Consequently the numerical 
Reynolds number R = LV/v ~ (8\Z2^) _1 ~ 90/(1000i/). 
This is very low and lower than most observed transitions 
between laminar and turbulent flows in the laboratory or 
terrestrial applications. 

The analytic solution in equation Q only holds if the 
fluid is Newtonian so that the shear stress can be described 
by the single number of the kinematic viscosity v. Because, 
Figure |A1| shows that the initial velocity profile changed 
strongly one also concludes that the fluid flow as modeled by 
SPH is that of a non-Newtonian fluid. So while it would have 
been interesting to think of SPH as solving effectively a the 
Navier Stokes rather than the Euler equations we see that 
this is not the case. The effective numerical viscosity is non- 
Newtonian and does not in general decrease with numerical 
resolution. 



3.4 Kelvin— Helmholtz Instability 

The Kelvin-Helmholtz instability occurs at the interface of 
two shearing fluids of different densities when velocity per- 
turbations perpendicular to the interface grow to eventually 
mix the layers. In the inviscid case this is understood an- 
alytically ( Chandrasckhar 1961) and the growth time scale 
for a sinusoidal mode of wavelength A between two fluids of 
density pi and p2 with a shear velocity v = t>2 — «i between 
them is 



2tt 



tkh 



(pi + P2)A 



(9) 



(PlP2) 1/2 v' 

This problem is typically setup with an adiabatic index of 
7 = 5/3 in a unit square domain x € {0, 1}, y G {0, 1} with 
(e.g. |Read et al.|2009") 



P,T,v x = 



pi,Ti,vi ly — 0.51 < 0.25 
P2,T 2 ,v 2 \y - 0.5| > 0.25 



(10) 



We choose p\ = 2, 5, 10 and p2 = 1 and the uniform pressure 
P2/7 = 3/5 which gives a sound speed of 1 in the low density 
surrounding medium. This standard setup then perturbs the 
interface with 



<5w y [sin(27r(x + A/2)/A) cxp(-(10(y - 0.25)) 2 
+ sin(27rs/A) exp(-(10(y - 0.75)) 2 )] 



(11) 



where we choose A = 1 and vary 5v y . So for our density 
contrasts the growth times are tkh ~ 2.12, 2.68 and 3.49 
for the density contrasts 1:2, 1:5, and 1:10, respectively. 

It would seem prudent to compare these results to the 
many investigations that recently have addressed the KH in- 
stability using SPH (Agertz et al.|2007||Wadsley et al.|2008 



Price|2008||Read et al.|2009||Hess and Springel|2009|> . How- 



ever, all of them chose a setup that Robertson et al. (20091 
showed to be ill defined. While all these studies where con- 
cerned with the issue of whether SPH can handle KH insta- 
bilities at all, they do not ask whether it actually converges 
to a correct solution. This can be seen e.g. in Read ct al. 
(2009) where their modified SPH solution compares already 
visually poorly to the corresponding Eulerian result. 
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Following |Robertson et ah| < |2009[ ) we opt for a more well 
defined setup for which they explicitly showed convergence. 
We modify the initial density, and velocity profile using the 
"ramp" function 



R(v) = 



1 



1 + exp[2(y - 0.25) /A„] 1 + exp[2(0.75 - y)/A y ] 



,(12) 



choosing A H = 0.05 so that p(y) — p2 + R(y)[pi — pi\ and the 
velocity shear is v x (y) — V2 + R(y)[vi — U2]. For the initial 
velocity perturbation we take v y (x) — 8v y sin(n7rx) setting 

0.01 and n = 2. 

Figure |3| compares results obtained with Enzo (Bryan 



SVy = 0.1 



an 



and Norman||1997 Bryan et al. 2001 O'Shea et al. 



2004), 



standard SPH and rpSPH for the two different initial veloc- 
ity perturbations. The improvement of rpSPH over SPH is 
dramatic. The billows grow at the correct rate and show 
a minimum of artificial small scale structure. The figure 
shows the highest resolution we have calculated. However, 
even with 120 2 particles neighbours one can obtain correct 
results using rpSPH. Also our choice of a high a is incon- 
sequential in rpSPH in this incompressible setup since the 
Balsara switch reduces it dramatically in practice. We left 
this high value just to show that one can get an accurate 
solution without having to adjust his viscosity parameter. 



3.5 Rayleigh Taylor instability instabilities 

Another classic test of a code's ability to handle subsonic 



perturbations is the RT instability (e.g Fryxell et al.||2000 



Stone and Gardiner|2007] |Stone et al.|2008[ ). SPH has been 
used to model supernova explosion previously and RT in- 



stabilities have been seen (Herant 



1994 1 as well as global 



convective instabilities ( Fryer|2004 1. The physical situation 
( Chandrasckhar 1961} here consists of a heavy fluid being 
supported by a lighter fluid against which initially are in 
pressure equilibrium with a constant acceleration (e.g. grav- 
ity). Remarkably, all idealised test cases that we are aware 
of use an initially unresolved contact discontinuity and con- 
sequently no converged results independent of the method of 
solution have so far been presented. Instead the differences 
between different reconstructions schemes, Riemann solvers, 
meshing, etc. all contribute to the final structures produced 
in these simulations. 

Similarly as the Kelvin-Helmholtz problem above, we 
chose here initial conditions that try to minimize the compu- 
tational requirements while yielding converged results. The 
two dimensional domain is setup with periodic boundary 
conditions along the x direction with x £ {0, 1/2} and re- 
flecting boundary conditions along y with y 6 {0.1,0.9}. 
To achieve the reflecting boundary conditions in Gadget- 
2 we set up the density distribution in 1/ £ {0, 1} and 
make all particles with y < 0.1 and y > 0.9 station- 
ary (-DSPH_BND_P ARTICLES compile option). These particles 
are not allowed to change their entropy or positions conse- 
quently they retain their initial pressure and density. Par- 
ticles that are trying to penetrate through these walls have 
their positions changed to be at the wall and y velocity vec- 
tors reversed. This is only a crude way of modeling reflect- 
ing boundaries with SPH but will suffice here to compare 
between SPH and rpSPH. 

In keeping with past literature we use an adiabatic index 
of 7 = 1.4, and set up the density at the top to be pi = 2 



t=3 
0.1 SPH 




t=4.5 
0.01 SPH 



t=4.5 
0.01 rpSPH 




Figure 4. Comparison of the final density distribution in a two 
dimensional Rayleigh Taylor test between standard SPH (left 
panels) and rpSPH (right panels) for two different initial veloc- 
ity perturbations Sv = 0.1 (top panels) and Sv = 0.01 (bottom 
panels). The unphysical "surface tension" of SPH prevents the 
growth of the instability entirely. rpSPH easily recovers the cor- 
rect behaviour. All simulations here used 500 2 particles and 70 
neighbours. 



and p2 — 1 at the bottom. So the density profile is p(y) = 
P2 + (pi - p2)/[l + exp(-2(y - 0.5)/A„)] with A y = 0.05 
in the cases presented here. The velocity perturbation is 
applied in y direction with v y (x,y) = 8v y (l + cos(87r(x + 
1/4)))(1 + cos(2/0.4 7r(y - l/2)))/4 and the y velocities are 
set to zero for y positions above 0.7 and below 0.3. The 
pressure is set to Po = P1/7 = 10/7 to give a sound speed of 
one at the interface and is set into hydrostatic equilibrium 
with the constant acceleration g — 1/2 in the negative y 
direction with P(y) = Po — gp(y)(y — 1/2). This gives a 
pressure difference of 60% between the top and the bottom 
of the domain. The smaller this pressure difference the more 
difficult it becomes for SPH to model it. Similarly, the initial 
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Figure 3. Comparison of the final density distribution in a two dimensional Kelvin-Helmholtz test. The velocity perturbation is a tenth 
of the sounds speed (0.1) in the top panels and one hundredth (0.01) in the bottom row. The left simulations show a grid based solution 
with 256 2 using the piecewise parabolic method. The middle row gives results for standard SPH and on the right are rpSPH calculations. 
This problem is set up care fully to not have any p hysical small scale structure and has been shown to converge with grid codes using a 
resolution exceeding 2500 2 l |Robertson et al.|2009| . rpSPH behaves significantly better in not breaking up from discretization noise. The 
bottom left shows how SPH fails to grow at the correct rate and is generally more noisy. All the particle based simulations used 1000 2 
particles and 100 neighbors with a = 3. 



velocity perturbation again should not be very much smaller 
than the sound speed in order to survive viscous damping 
before a growth time of the instability. 

We present the results for this test for velocity pertur- 
bations of Sv — 0.1 and Sv = 0.01 in Figure [4] 

The differences are dramatic. Where SPH fails com- 
pletely to see growth of the instability rpSPH gives the ex- 
pected behaviour for both perturbation strengths. 

That rpSPH is dramatic improvement over Morris' for- 
mulation despite only differing in one index is seen in Fig- 
ure [5] There we give a Rayleigh Taylor problem at low res- 
olution of 100s50 particles and a density ratio of 10 as fur- 
ther discussed in section |4~T] All parameters were the same. 
A courant factor of 0.2 is used, a neighbour number of 40, 
a = 1.5, Balsara switch is on, and the initial velocity per- 
turbation amplitude is 0.1. Clearly our formulation is more 
stable lending support to our discussion on the different er- 
ror properties of the two discretisations given above. 



3.6 Shock tubes 

3.6.1 Sod shock tube 

So far we have tested our new formalism only in very weakly 
compressible situations. We will use the classic Sod shock 



tube (|Sod||1978[ ) to compare rpSPH to standard SPH here. 
Rosswogj ( |2009[ ) recently, gave the results for varying viscos- 
ity prescriptions and including artificial conduction terms. 
We change the setup only slightly. The left state has a den- 
sity and pressure of unity while the right state has a quarter 
of the density and a pressure of 0.1795. This test is evolved 
with an adiabatic index of 7 = 1.4 and we set it up as a 
two dimensional problem with equal mass particles in a box 
that extends from zero to ten in x and zero to one in the 
y-direction. We choose 40 rows of particles in the y direction 
and vary the spacing along x to achieve the given densities 
using a total of 200 2 particles which are initially at rest. Ad- 
ditionally we set the interface to be at x = 3 and smooth 
it with an exponential ramp such that all hydrodynamic 
variables are given by r + (1 + exp(2 * (x — 3)/Sx))~ 1 (l — r) 
where we take Sx — 0.05 and I and r denote the left and right 
states. We employ periodic boundary conditions which gives 
us another interface at x = 10 which has the reversed left 
and right states but has an initially discontinuous state. This 
will give us the opportunity to show the difference between 
smoothed interfaces and artificially sharp to be visible in one 
figure. For the first results we present we have used 80 neigh- 
bors and an artificial viscosity parameter of a = 3 for both 
the SPH and the rpSPH calculation. Both employ the Bal- 
sara switch to limit the viscosity which will play no role here 
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Figure 6. Clockwise the density, rr-velocity, pressure and specific internal energy are given for the classic Sod shock tube jSod|1978| l at 
time t = 1. The simulations used 200 2 particles, 80 neighbours with a tolerance of one, a = 3 for both the give the rpSPH solution (filled 
diamonds) and standard SPH (open circles). We plot every 50th particle of the 40,000 employed. In the pressure panel (bottom right) 
we also show a references solution computed with a Eulerian code using the piecewise parabolic method. The agreement is quite good. 
rpSPH handles the contact discontinuities much better and the "blip" seen in SPH is absent in the one form the smoothed interface 
(x ~ 3.8) while it is much reduced in the one that was initially sharp (x ~ 9.3). 



because there is no rotational component to the flow. Fig- 
ure [6] summarizes the hydrodynamic state variables at time 
t = 1. To first approximation one gets identical results with 
the new formalism as compared to the standard approach. 
Linear momentum is not conserved in rpSPH and we find 
a linear excess velocity of (—105, —3 x 10 -5 ) so per particle 
an error on the velocity of 0.003 in the x-direction and a 
completely negligible component along the y-direction. This 
is at a time when the r.m.s. velocity is ~ 0.46 so just slightly 
above one half of a per cent error in the dominant a;- velocity. 
SPH has poor behaviour at the contact discontinuities. For 
both the one originating from the initially smoothed and the 
the discontinuous interface at the right boundary. Both con- 
tacts at x ~ 3.8 and x ~ 9.3 are better captured by rpSPH. 
The Sod shock tube has few features and it is reassuring 
that using as many as 200 2 particles can give an excellent 
answer. There are only slight differences in how rpSPH han- 
dles one dimensional shock tubes. We will discuss one very 
popular application taken from a cosmological context after 
testing a very strong shock next. 



3.6.2 Strong Shock 

Here we give another test of a much stronger shock than 
the one by sod. This one has a Mach number close to one 
hundred. We also use the chance to compare this to the dif- 



ference formulation studied by Morris (1996a). The density 
and pressure are (1, 6.6 x 10 4 ) on the left and (1/5, 1) on the 
right. This is very similar to the one studied by |Pfrommer| 
et al. ( 2006 I and is well known to work well with standard 
SPH. Here we use 35 neighbours, a = 4 and 5000 particles. 
This is a good example where one can make rpSPH and the 
Morris formulation give unphysical results. These methods 
require the pressure gradient to be resolved. So if you start 
with completely discontinuous left right states one will get 
unphysical waves giving unexpected results. However, this is 
not a shortcoming of the method but simply are errors that 
come from not resolving the initial conditions. We again use 
the ramp function from above with a width of 4 in this very 
long domain ranging from to 500 in x and to 10 in y. 

We cannot confirm Morris' claim that his formulation 
gives large post shock oscillations in this method and suspect 
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Figure 7. Clockwise the density, x-velocity, pressure and specific internal energy are given for a strong shock with Mach number ~ 100 
at time t = 0.5. We used The simulations used 5000 particles arranged as ten rows in y in an elongated domain 0, 500, 0, 10, 35 neighbours 
with a tolerance of one, a = 4. We give the result obtained with Morris' formulation (open circles) and rpSPH. We plot every 50th 
particle of the 5,000 employed. Both approaches handles the contact discontinuities much better and the "blip" seen in standard SPH is 
gone. rpSPH gets the correct density jump of a factor of 4 better than the Morris formulation. 



that he may have set up discontinuous initial conditions. We 
can see that our new formulation performs somewhat better 
than the Morris formulation as it does not overshoot the 
analytical density jump of 4 raising the density from 0.2 to 
0.8 in Figure [7] Otherwise both approaches work fine and 
have no problem in modeling strong shocks and evolving it 
for large distances. 



3.7 Sedov-Taylor Blast Wave 

Another particularly strong shock is formed in the Sedov- 



Taylor blast wave (Landau and Lifshitz 19591 presenting 



a difficult test problem for incompressible hydrodynamics 
codes. One the one hand it is a self similar solution which 
makes it insensitive to how exactly one sets it up as long as 
one evolves the system for a very long time. On the other 
hand it is the solution for a point explosion. For a given finite 
resolution, however, there is no unique way of specifying the 
initial conditions. Here is where exact momentum and en- 
ergy conservation is very helpful as one can set up the initial 
conditions at will and even if one were to make very large 
errors in the time evolution the method will still arrive at 
the self similar solution. In conservative grid codes this still 



can lead to aspherical solutions if one did not resolve the 
spherical central hot region. Since SPH however uses spheri- 
cal kernels one can get away sometimes even by just heating 
one single particle (Springel and Hernquist 2002]) . This is 
very useful in applications such as galaxy formation simu- 
lations where one is always far from resolving the relevant 
length scales of an explosion. On the other hand any physics 
that were to occur at a scale of the shell thickness would be 
impossible to resolve in such a single particle energy ejection. 
For rpSPH and the Morris formulation we need to resolve 
the pressure gradients in the initial conditions as we saw 
above in the strong shock setup. 

We set up a square lattice of particles with 300 particles 
on a side in the unit domain. For resolved initial conditions 
we set a spherical region in the center of radius r = 0.1 with 
the same ramp function as above using a width of 0.1, to 
have a sound speed of one for an adiabatic index of 7 = 5/3. 
For both simulations we used a Courant number of 0.2 (0.1 
in Gadget), 80 neighbors, artificial viscosity a — 2.5 and had 
the Balsara switch off. Figure [8] shows that there potentially 
is also an advantage to rpSPH simulations when modeling 
shocks. We have failed to get standard SPH to give a stable 
correct density jump of (7 + l)/(7 — 1) = 4 for our setup. 
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Also the three dimensional versions shown by |Springel and| 
Hernquist ( 2002 1 always seem to be too low by as much as a 



Figure 5. Density in a Rayleigh-Taylor test with a ten times 
heavier fluid on top than the bottom. Left is the Morris formula- 
tion and on the right is rpSPH. Both Morris SPH and rpSPH runs 
used an initial uniform grid of 50 by 100 particles with varying 
masses to describe the higher density for the interface and top 
fluid. The Morris formulation is unstable in this problem while 
rpSPH behaves as expected. 
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Figure 8. Comparison of the shock profile in the density between 
SPH and rpSPH simulations. The analytic solutions demands a 
density jump of four which rpSPH and the Morris formulation 
stay close to over the course of the calculation. Standard SPH 
falls a little short, however, it is insensitive to choices of initial 
conditions and is a better choice if one cannot afford sufficient 
resolution to model the exploding region. 



factor of two. However, standard SPH is much less sensitive 
to how one sets up the initial conditions and performs much 
better at low resolutions. 

These strong shock problems can work but clearly are 
not the biggest strengths of rpSPH. However, at this point 
we simply reused the artificial viscosity prescription which 
was designed for standard SPH. We believe it is likely that 
one can find an alternative formulation for the artificial vis- 
cosity that fits better into the rpSPH discretisation which 
may improve its behaviour for highly supersonic conditions. 
Until a better artificial viscosity prescription is designed one 
may opt to switch between standard SPH and rpSPH based 
on the local divergence. We have successfully applied this 
strategy by using a switch that evaluates the standard SPH 
sum if —hi -i-Vi> 3c Si i and the rpSPH sum for less strongly 
convergent flow. Here hi and c Sl i denote the smoothing 
length and the current sound speed of the particle. This 
formulation is robust in all our tests. 



3.8 Cosmological Integration of the Santa 
Barbara Galaxy Cluster 

In 1995 a comparison project was initiated that aimed to 
compare all numerical cosmology codes at that time for rel- 
evant realistic initial conditions. The study focused on three 
dimensional calculations of the formation of a galaxy clus- 
ter in the standard CDM scenario of structure formation. 
The choice was a setup which does not include any other 
physics than cosmological hydrodynamics with an ideal gas 
equation of state (often referred to as adiabatic simulations 
despite the entropy generation in shocks). The study pro- 
duced a detailed report in ( Frenk et al.[1 999 F99, herafter). 
One of the most surprising findings of the study was that 
while there was very good agreement between the six dif- 
ferent SPH codes used in the study they did not agree with 
the solutions of the grid based codes. The central entropy of 
the simulated galaxy cluster differed markedly between the 
grid and SPH codes. In particular, the only AMR code in 
the study by Bryan and Norman ( 1997 1 which is now called 
Enzo ( |Bryan et al.|2001||O'Shea et al.|2004[ ) found a flat en- 
tropy core while the particle codes found the central entropy 
to continue to rise towards smaller radii. Note that there has 
been a significant debate on what the correct solution may 
be and potential sources between the differences between 



grid and particle based methods (|Springel|2005l IDolag et al. 



2005l | Kawata et al.|2 009 Wa dsley et al. 2008, Mi tchell et al 
2009| |Agertz et al.|2007| Springel|2010[ ) is the real reason for 
the difference. 

We do not attempt a resolution study here but sim- 
ply show how an rpSPH solution compares to an SPH run 
with otherwise identical parameters and the solution derived 
with a cosmological AMR code. We use 128 3 gas as well as 
dark matter particles for the particles based approach. For 
the AMR code we again use Enzo already used in F99 using 
128 3 dark matter particles a root grid of 128 3 cells and seven 
additional refinement levels. Refinement is based on density 
thresholds in the baryons and dark matter component. The 
viscosity parameter in the particle based runs was a — 3, 
and a neighbor number of 300 with tolerance of one was 
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Figure 9. Logarithm of temperature in a two dimensional slice at redshift zero derived form the initial conditions of the Santa Barbara 
Galaxy Cluster comparison project ( |Frenk et al.|1999) . The plots are 64 Mpc across and show 10 contours within the interval from 30,000 
to one hundred million degrees Kelvin. The SPH and rpSPH simulations used 128 3 particles and 300 neighbors with a tolerance of one, 
a = 3. The agreement is excellent. Subtle differences in the postshock temperature can be seen. rpSPH seems clearly robust enough for 
large scale cosmological integrations. 



used. The initial redshift is z = 63 and the initial tempera- 
ture 10 9 Kelvin. 

A two dimensional slice through the temperature field is 
given in Figure [9] comparing SPH to rpSPH. Only relatively 
subtle differences are found. There are perhaps slightly more 
small scale features visible in the rpSPH calculation and 
some slight differences in the post shock gas in the main 
cluster are visible. Slightly more shocking occurs at larger 
radii towards low density voids in the rpSPH vs. SPH cal- 
culation. 

Figure [To] compares the solutions using spherically aver- 
aged radial profiles as described in F99. The entropy profiles 
of rpSPH agree better with the AMR than the SPH results. 
The rpSPH solutions shows the lowest central densities of 
the three methods and agrees better with the slightly shal- 
lower density profile of the grid code. Clearly the differences 
between all three, however, are rather subtle given that one 
evolved this system for 13 billion years which are many tens 
of sound crossing times for the central part of the resulting 
cluster. 

Both SPH and rpSPH simulations have a final linear 
momentum corresponding to 4.1 and 3.2 km per second per 
gas particle, respectively. The difference vector between the 
final gas momenta of the simulation has a magnitude of 2.6 
km/s per gas particle. This is a difference of order one half 
of a percent of the mass weighted mean r.m.s. velocity of 
~ 400 km/s. Obviously, giving up the strict linear momen- 
tum conservation in our equation of motion has not lead to 
any noticeable difference in this measure but has improved 



the comparison with results from adaptive mesh refinement 
codes. 

However, this particular application is relatively easy 
as dark matter dominates the gravitational potential. As we 
will show further below rpSPH is quite easy to break with 
self-gravitating fluids. 



4 FURTHER BENEFITS OF RPSPH 
4.1 Variable masses 

Next we demonstrate that using our pressure force discreti- 
sation give another very important advantage. Simulations 
with drastically varying particle masses continue to give cor- 
rect results. This is markedly different compared to previ- 
ous SPH simulations employing particle splitting. The latter 
only worked reasonably as long as different particle masses 
were very well separated spatially. As an explicit example we 
revisit the Rayleigh-Taylor problem from above. This time 
we initialize particles on a uniform lattice and model the 
density contrast by changing the particle masses according 
to the density profile. We employ a density at the top ten 
times the one of the one at the bottom fluid to demonstrate 
that this is not just a marginally better aspect of rpSPH. 
Figure [TT] compares the SPH and the rpSPH solution again 
for 5v = 0.01. Instead of showing the density field we show 
the particles painted by circles and colored by their density. 
This gives us an opportunity to see that rpSPH does not 
show any clumping instability while it is severe for SPH. 
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Figure 10. Spherically averaged profiles of density, temperature 
and entropy of the Santa Barbara Galaxy Cluster comparison 
project | |Frenk et al.|1999] l The SPH and rpSPH simulations used 
128 3 particles and 300 neighbors with a tolerance of one, a = 
3, and are compared to results from an AMR calculation using 
the piecewise parabolic method (enzo-PPM). The agreement is 
excellent. Its interesting that rpSPH gives higher central entropy 
as standard SPH more comparable to the grid code as this has 
been the subject of much discussion in recent years. Note also 
that the entropy profiles of rpSPH and the AMR calculation agree 
better well away from the central resolution limited part at radii 
between 0.1 and 1. 




Figure 11. Density in a Rayleigh-Taylor test with a ten times 
heavier fluid on top than the bottom. Both SPH and rpSPH runs 
used an initial uniform grid of 100 by 200 particles with varying 
masses to describe the higher density for the interface and top 
fluid. The clumping instability of standard SPH creates the pat- 
tern on the left while rpSPH works perfectly fine. Being able to do 
accurate calculations even with largely disparate particle masses 
at low resolution will be an enormous benefit in many situations 
where large density contrasts are observed but one needs to retain 
high accuracy in the low density regions. 



4.2 Formulation for Magnetic Forces 

There have been many attempts to implement ideal MHD 



into the standard SPH formalism (e.g. Gingold and Mon- 
aghan] |T977| [Phillips and Monaghan 1985 ; Price and Mon- 
aghan|[2004| |Dolag and Stasyszyn||2009| ) In our experience 
rpSPH performs in the hydro part better than the Morris 
formulation. The latter has been used in implementing ideal 
MHD into SPH (|Morris||1996b| |P~rice and MonaghanpXM] 
Dolag and Stasyszyn||2009 |. Therefore, we expect that our 
new discretization may be of use in this case as well. 

A symmetric conservative form of the Lorenz force 
is generally implemented using the magnetic stress tensor 
( |Phillips and Monaghan|1985[ |, 

1, 



\B l \ 2 S kl 



(13) 



In this run we only use 100 by 200 particles demonstrating 
that rpSPH handles the Rayleigh-Taylor problem very well 
at small initial perturbations and low particle numbers. 

The Sedov- Taylor blast wave above was also carried out 
with a staggered grid of particles of varying mass and re- 
tains a nice spherical shape despite the multiple squares 
introduced in the staggered "mesh" of the initial particle 
distribution. 

It is of great interest for a method to be stable under 
largely varying particle masses. If it is one can use particle 
splitting likely without worrying too much of how to place 
the new particles and keep them separate from particles with 
different masses (e.g. Kitsionas and Whitworth||2002 1. 



So the acceleration from the magnetic fields on the i-th par- 
ticle is then written as 



dv, 



JV 

1 x - 



/to 



Mi - Mi 



P7 



In the limit where we only have forces from the magnetic 
pressure gradient only the diagonal of the tensor has terms 
which are Bf/2 and we recognize this discretisation as ex- 
actly the form of equation |2| above. So also the Lorentz 
force lends itself to be discretised following our new ap- 
proach. We split the force into the tension component and 
the magnetic pressure component without loss of generality 
into 
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= -V 



B 2 (B • V)B 



871 



+ 



4tt 



(14) 



In both terms one takes spatial derivatives and hence is al- 
lowed to subtract a constant. Choosing Bi ■ Bi for the first 
one and Bi for the second one avoids finite pairwise forces 
between particles in regions of constant field. So the isotropic 
part becomes 



dvi 

~dT 



N 



3 = 1 



£?2 p2 



vw, 



to which the tension force is added 



+ -E 



Bi ■ (Bj — Bi 



VW. 



(15) 



(16) 



3=1 



Both the terms simply add to the accelerations in the 
force calculation. All that is left to do is to replace the fastest 



signal velocity with the one given in equation (46) of Price 
land Monaghan ( 2004 1 and one has a MHD implementation 



of SPH. On some simple initial tests this formulation seems 
to work quite well even without any regularization technique 



(e.g. Dolag and Stasyszyn 2009 and references therein) or 



artificial B field dissipation. A full exploration of the perfor- 
mance of this discretization though is beyond the scope of 
this contribution. 



5 HOW TO BREAK RPSPH 

From the tests above it is clear that at least in the weakly 
compressible limit rpSPH is a very useful improvement over 
the standard formulation. However, giving up momentum 
and energy conservation is a big price to pay for those ad- 
vances and rpSPH cannot possibly replace standard SPH 
in all problems of interest. rpSPH should be easy to break 
at low resolutions when one does not resolve the pressure 
gradients adequately. We now give an illustrative example 
that makes rpSPH give very bad results which shall serve 
as a cautionary note and things to look for when applying 
rpSPH. 

The Evrard collapse of a cold gas sphere ( Evrard|1988 1 
has been extensively used for verification of astrophysical 
SPH codes (e.g. |Wadsley et al.||2004| |Springel||2005[ ). In the 
version that is part of the Gadget distribution it is realized 
with only 1472 equal mass particles. A centrally concen- 
trated cloud of cold gas collapses, bounces and eventually 
virializes. Vacuum boundaries are assumed. Very clearly this 
is only meant to investigate energy conservation of the code 
and is a simple test running in seconds on ones laptop. 

Figure [12] shows the standard SPH solution together 
with a completely failed solution of rpSPH. We should not 
that as we increase the resolution rpSPH does converge to 
the same solution as standard SPH. However, this is a clear 
example where too little resolution coupled with a solver 
that does not conserve momentum or energy will fail com- 
pletely. 

Figure [12] also gives the energies in the Evrard collapse 
for which we formally added a zero and label it rpSPH+0. 
The term added to the momentum equations right hand side 
is P/pVl — 0, which in SPH reads (e.g. Ritchie and Thomas 
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Figure 12. Energies in the Evrard gas collapse test. 
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This term which tends to zero when the density and pressure 
gradients are well resolved. For our poor resolution setup, 
however, we see that it gives better energy conservation. 
We do not recommend this formalism over rpSPH, however, 
since it still has the problems with contact discontinuities 
and weakly compressible flows. However, it highlights that 
as for any numerical study resolution studies are crucial and 
that rpSPH will likely be of little use when one cannot af- 
ford sufficient resolution for the particular problem one is 
interested in. 



6 CONCLUSIONS 

We have presented a novel discretization of the pressure 
equation for the smoothed particle hydrodynamics, which 
we call rpSPH that removes the local pressure from the 
scheme and only considers pressure gradients. This method- 
ology avoids the clumping and banding instability, artifi- 
cial surface tension, unphysical particle noise, dramatically 
reduces inherent shear viscosity and numerical dissipation, 
and allows to realistically evolve density distributions sam- 
pled even with disparate particle masses. We have discussed 
a large number of test in all of which our new discretiza- 
tion outperforms the traditional SPH results. While our ap- 
proach is not manifestly momentum conserving and easy 
to break wiht self-gravity and or low resolutions it clearly 
seems much more accurate than previous approaches to La- 
grangian hydrodynamics using SPH. Since our formulation 
is more accurate and requires as little as one line of code to 
be changed in previous implementations we do believe it to 
likely to be useful. The caveat of rpSPH remains that if one 
knows that one cannot afford to resolve the pressure gradient 
in ones initial conditions that because it is not momentum 
conserving it can give very wrong results. Fortunately, a res- 
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olution and convergence study can reveal whether one is in 
this limit. 

In summary, some of the biggest shortcomings of SPH 
can in some circumstances can be overcome if one gives up 
the idea of applying equal but opposite forces to particle 
pairs. While the latter is what happens physically to the 
atoms or molecules making up the fluid it is simply incor- 
rect for Lagrangian fluid elements the particles are meant 
to represent. Physically it also does not make sense to in- 
troduce repulsive forces for two spatially separated points 
even when there is no pressure gradient between them. To 
require such symmetry between particles neglects that they 
are spatially separated and that the gradient of the pres- 
sure field is different at the two locations in general. Our 
new discretization avoids these unphysical forces and allows 
the SPH particles to behave as Lagrangian volume elements 
recovering fluid behaviour in a large number of tests. 

We have successfully used a fifth order spline kernel giv- 
ing smaller errors on the uniform shear problem and the 
Rayleigh- Taylor problems discussed above. Consequently, 
we believe that further improvements to rpPSPH should be 
possible in the future. 

We have also studied multiple forms of discretising the 
specific internal energy equation 



The simplest version that we successfully applied to some of 
our test problems is given by 

3=1 

While we prefer the entropy formulation this form here may 
be useful for codes that start from an internal energy for- 
mulation. 

While standard SPH is conservative it fails to correctly 
capture fluid instabilities and shows large non-Newtonian 
viscosity. rpSPH, on the other hand is more accurate, but 
is not inherently momentum or energy conserving. Conse- 
quently it is a useful modification to the SPH algorithm 
when one is studying problems where one can afford to re- 
solve the relevant pressure gradients and the density field. 
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APPENDIX A: SHEARING FLOWS OF 
UNIFORM DENSITY WITH STANDARD SPH 

While rpSPH overcomes the problems found for shearing 
flows this appendix is more of historical interest. However, 
some readers may find it worthwhile since to our knowl- 
edge the surprisingly large shear viscosity and its erratic 
behaviour with increasing particle numbers has not been 
documented previously. 

A simple two dimensional setup uses an adiabatic index 
of 7 = f .4 in a unit square domain x £ {0, f}, y G {0, f} 
with periodic boundary conditions. All particles are set up 
on an exactly square lattice with a density of p(x, y) = f so 
that the initial density estimate from the SPH kernel in fact 
gives a density estimate of unity to better than four parts in 
one thousand. We then add different velocity perturbations 
to this uniform distribution. We set the pressure to Po = p/j 
to have a sound speed of unity. For the first tests here we 
only use 50 2 particles as there are no features to resolve. In 
all cases we evolve to time t = 4. 

First we start with no velocity perturbation. I.e. a com- 
pletely static uniform density distribution evolved over four 
sound crossing times. Using 30 neighbors the density esti- 
mate by all particles is 1.00345. For different neighbor num- 
bers this fluctuates around I and is close enough. We will 
use 30 neighbors for most of the rest of this section. Ini- 
tially all velocities are zero yet after four crossing times 
we have an r.m.s. velocity v rms = U v % + v v ~ 0.01 
with no obvious preferred direction. This results is obtained 
for the typical viscosity value of a = 1. For lower values 
this random noise increases to v rms ~ 0.036 for a = 1/10. 
So clearly even under the most quiet conditions imagin- 
able, a uniform density in pressure equilibrium, we could 
not represent velocities of order a few percent of the sound 
speed. Now let us perturb the velocity along the x direc- 
tion and set it to a uniform value of 1. This should be 
exactly identical to the previous setup given that SPH is 
formulate d to be Galilean invariant . Now the random noise 
vrms = ^/iV 'EiVpO* - I) 2 + v y ~ 0.009 for a = 1 and 

again v rms « 0.035 for a = 1/10. 

Now for our next experiment with this uniform density 
setup we use v x (y) = 5v y cos(27ry) with 5v y = 1/2. This 
shear flow setup gives an average kinetic energy of 1/8. Af- 
ter only 4 sound crossing times (or two crossing times of the 
fastest particles) the mean kinetic energy of particles has 
decreased by 15% and the r.m.s. velocity in y direction is al- 
ready » 0.063 when using a viscosity parameter of a = 1/10. 
Using the standard value a = 1 we have a lower r.m.s. ve- 
locity in the y direction of w 0.032 yet at the same time the 
total kinetic energy has decreased by a as much as 27 percent 
suggesting that the standard value does convert unaccept- 
able levels of the shear into heat. This hardly is inviscid flow! 
Again for a — 1 but 200 2 particles which allow the shear to 
be better resolved one would hope for less dissipation. Yet 
we find that the total kinetic energy still decreases by 30% 
and the r.m.s. velocity in the y direction becomes « 0.021. 
We have also run this test with 300 2 particles and find that 
the kinetic energy dissipation is approximately independent 
of resolution up to this particle number. The kinetic energy 
lost after two crossing times was 30.3% and the final r.m.s. 
velocity fluctuations in the y-direction was w 0.019 i.e. a 
fiftieth of the sounds speed. The latter velocity dispersion 
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became as high as a thirtieth during the first time inter- 
val t ~ 0.2, decreased and then stayed stable afterwards at 
« 0.02. The maximal vertical velocities are as much as one 
tenth of the sound speed for a problem which should not 
develop any perpendicular velocities. 

Perhaps using a lower viscosity parameter could help? 
So with a — 0.1 and 200 2 particles this indeed gives lower 
overall dissipation of the kinetic energy of "only" 16% but 
then gives random motions perpendicular to the shear of 
« 0.030. Figure [2] summarizes the loss of kinetic energies 
for simulations with 50 2 and 100 2 particles for 30 and 48 
neighbors and the artificial viscosity parameters of a — 0.1 
and a = 1. 

From the figure it is very clear that using more particles 
actually leads to more dissipation. How is one supposed to 
carry out a resolution study when the numerical dissipation 
can keep increasing when using more and more computa- 
tional resources? 

A histogram over all particles showing their x velocity 
as a function of their y coordinate in Figure 3 reveals how 
strongly the shear viscosity turned the initial sinusoidal per- 
turbation into flattened extrema with linear profiles between 
them. This graph looks similar for different viscosity values. 
The bottom panel of Figure 3 visualizes the particles making 
the ones that received the entropy clearly visible. The lowest 
(initial) value one should expect is white and in fact below 
the minimum on that image Pjp 1 = I/7 ~ 0.7143. How- 
ever, one can see clear bands in the places where one finds 
the largest gradients in the shear velocities. The clumping 
instability is clearly visible through the bunching of entropy 
values in the plot. 



It is worth noting that the Balsara switch (Balsara 



1995) which is designed to limit this shear dissipation in- 



deed helps. Without it we find that after two crossing times 
48% of the kinetic energy are already artificially dissipated 
in the 200 2 test with the sinusoidal shear at Mach one half 
and a uniform density. These 48% are to be compared to the 
30% which were dissipated using the Balsara switch. 

The viscosity limiter implemented in Gadget-2 does not 
influence the results here. Also decreasing the Courant factor 
by an order of magnitude can change the exact amount of 
dissipation but does in general not decrease it appreciably. 

That the effective shear viscosity changed little with 
increasing particle number is very unfortunate. However, if 
we increase the neighbor number employed to 60 neighbors 
the effective shear viscosity drops dramatically and only 1% 
of the total kinetic energy is artificially dissipated over the 
same time interval. However this comes at the price of par- 
ticles clumping into bands through the well known tensile 
instability. As discussed in some detail by ( [Read et al.|2009[ ) 
the amount of clumping is specific to the kernel choice. 

One compromise for this uniform shear problem is a 
neighbor number of 48 which leads to banding and only 
about 3% kinetic energy dissipation in the two crossing 
times. Simple scaling implies then a choice of 48 3 ^ 2 ~ 333 
neighbors for three dimensional calculations. A number 
much larger than typically employed. 

We have evolved this same test to many more crossing 
times and find that at the larger neighbor numbers dissipa- 
tion simply occurs later but in fact looks qualitatively just 
like in the low neighbor number case. This enables now a 
further discussion of the origin of this artificial shear viscos- 




0.74 



0.73 



0.72 



Figure Al. Uniform density shear test. Top: The solid line gives 
the amplitude of the uniform initial x velocity modulated as a 
function of y. The underlying histogram is the particle distribu- 
tion after four sound crossing times when using a viscosity param- 
eter a = 1. It looks slighlty better for the lower artificial viscosity 
parameter of a = 1/10 but leads to larger vertical velocity pertur- 
bations. Bottom: Entropy of the particles in a 200 2 run showing 
bands of the material that received the entropy (P/p 7 ) in the 
200 2 run with a = 1. The values are plotted as squares with side 
lengths of half the SPH smoothing length. 



ity. In the top panel of Figure |A1| we can see how the ex- 
trema in the velocity are clipped by the viscosity and that 
the particles positions which were one a regular lattice in 
the y direction spread out and led to banding. This in large 
part comes from the non-zero y velocity the particles obtain 
leading them to artificially mix into regions perpendicular to 
the velocities they have. So small fluctuations in the pressure 
forces enable a coupling taking energy from the x velocity 
to stimulate motions in the y direction. Particles then artifi- 
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cially mix into regions of the flow where they start to inter- 
act with fluid parcels of different shear velocities triggering 
the artificial viscosity which then tries to damp this noise. 
This is why smaller artificial viscosity parameters will lead 
to larger r.m.s. y velocities. This also explains why increas- 
ing the neighbor number delays this artificial shear viscosity 
in that it decreases the amplitude of the forces leading to 
the y velocities. 

Interestingly the uniform density tests presented here 



are similar as the one presented by Monaghan ( 2006 1 in 



terms of the velocity profile and in the sense that it is a 



low mach number flow. |Monaghan (20061, however, chose 
to pick 7 = 7 making the EOS very stiff. While this may 
be a useful trick to model incompressible flow with SPH it 
is not something we repeat here since for most applications 
in astrophysics we have 1 ,$ 7 5/3. However, even with 
the stiff equation of state his Figure 2 shows the same clip- 
ping of the maximal velocity amplitudes as our Figure |A1| 
after only one crossing time. This agrees with our findings 
that only for short time scales (as compared to the crossing 
time) the shear viscosity may be negligible. We just differ in 
the interpretation of whether this is an acceptable level of 
dissipation or not. 



APPENDIX B: MODIFYING GADGET-2.0.4 TO 
RPSPH 

For the convenience of other researchers we give the details 
of what to do to convert Gadget-2.0.^] to take advantage 
of the rpSPH discretization. In hydra, c find the line that 
reads 

hfc = hfc_visc + P [ j ] . Mass * ( p _o ver _r ho 2 _i * dwk.i 
+ p _o vcr _r ho 2 _j * dwk_j ) / r ; 

and change it to 

hfc = hfc_visc+P[j ]. Mass/SphP[j ]. Density* 
(SphP [j ] . Pressure— pressure) / SphP [j ] . Density* 
dwk.i / r ; 

and the conversion is complete. 

Another form that fits more closely to the artificial vis- 
cosity prescription useful for problems with large density 
gradients is 

hfc = hfc_visc+P[j ]. Mass/SphP[j ]. Density* 
(SphP [j ] . Pressure — pressure)/ SphP [j ] . Density* 
(dwk_i+dwk_j )/2/r ; 

In order to keep standard SPH for very strong shocks 
matching the standard viscosity implementation one may 
choose to keep both lines but preface the latter with 

if (— h_i*divVel < 3 . * soundspeed _i ) 
, or other criteria that trigger at strong shocks. 

This paper has been typeset from a TgX/ F/TjtX file prepared 
by the author. 



1 http: //www.mpa-garching .mpg. de/gadget/ 
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